﻿using System;
using System.Text;
using System.Drawing;
using System.Buffers;
using System.Collections;
using System.Collections.Generic;
using System.Runtime.InteropServices;

public static partial class NativeAOT
{
    [UnmanagedCallersOnly(EntryPoint = "tlvs")]
    public static unsafe int tlvs(IntPtr t_ptr, int n, IntPtr b_ptr, IntPtr x_ptr)
    {
        double* t = (double*)t_ptr.ToPointer();
        double* b = (double*)b_ptr.ToPointer();
        double* x = (double*)x_ptr.ToPointer();

        return tlvs(t, n, b, x);
    }

    /// <summary>
    /// Toeplitz方程组
    /// </summary>
    /// <param name="t">t[n]存放n阶T型矩阵中的n个元素。</param>
    /// <param name="n">方程组阶数。</param>
    /// <param name="b">b[n]存放方程组右端的常数向量。</param>
    /// <param name="x">x[n]返回方程组的解向量。</param>
    /// <returns>若=0则表示失败；若>0则表示正常。</returns>
    public static unsafe int tlvs(double* t, int n, double* b, double* x)
    {
        int i, j, k;
        double a, beta, q, c, h;
        double* s = stackalloc double[n];
        double* y = stackalloc double[n];

        a = t[0];
        //if (Math.Abs(a) + 1.0 == 1.0)
        if (EQ0(a))
        {
            return (0);
        }
        y[0] = 1.0;
        x[0] = b[0] / a;
        for (k = 1; k <= n - 1; k++)
        {
            beta = 0.0; q = 0.0;
            for (j = 0; j <= k - 1; j++)
            {
                beta = beta + y[j] * t[j + 1];
                q = q + x[j] * t[k - j];
            }
            //if (Math.Abs(a) + 1.0 == 1.0)
            if (EQ0(a))
            {
                return (0);
            }
            c = -beta / a;
            s[0] = c * y[k - 1];
            y[k] = y[k - 1];
            if (k != 1)
            {
                for (i = 1; i <= k - 1; i++)
                {
                    s[i] = y[i - 1] + c * y[k - i - 1];
                }
            }
            a = a + c * beta;
            //if (Math.Abs(a) + 1.0 == 1.0)
            if (EQ0(a))
            {
                return (0);
            }
            h = (b[k] - q) / a;
            for (i = 0; i <= k - 1; i++)
            {
                x[i] = x[i] + h * s[i];
                y[i] = s[i];
            }
            x[k] = h * y[k];
        }
        return (1);
    }

    /*
      int main()
      { 
          int i;
          double x[6];
          double t[6]={6.0,5.0,4.0,3.0,2.0,1.0};
          double b[6]={11.0,9.0,9.0,9.0,13.0,17.0};
          if (tlvs(t,6,b,x)>0)
          for (i=0; i<=5; i++)
              cout <<"x(" <<i <<") = " <<x[i] <<endl;
          return 0;
      }
    */
}
